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ABSTRACT 

Evolution of stellar bars in disk galaxies is accompanied by dynamical instabilities and secular 
changes. Following the vertical buckling instability, the bars are known to weaken dramatically and 
develop a pronounced boxy/peanut shape when observed edge-on. Using high- resolution iV-body 
simulations of stellar disks embedded in live axisymmetric dark matter halos, we have investigated 
the long-term changes in the bar morphology, specifically the evolution of the bar size, its vertical 
structure and exchange of angular momentum. We find that following the initial buckling, the bar 
resumes its growth from deep inside the corotation radius and follows the Ultra-Harmonic resonance 
thereafter. We also find that this secular bar growth triggers a spectacular secondary vertical buckling 
instability which leads to the appearance of characteristic boxy/peanut/X-shaped bulges. The secular 
bar growth is crucial for the recurrent buckling to develop. Furthermore, the secondary buckling is 
milder, persists over a substantial period of time, ~ 3 Gyr, and can have observational counterparts. 
Overall, the stellar bars show recurrent behavior in their properties and evolve by increasing their 
linear and vertical extents, both dynamically and secularly. We also demonstrate explicitly that the 
prolonged growth of the bar is mediated by continuous angular momentum transfer from the disk 
to the surrounding halo, and that this angular momentum redistribution is resonant in nature — a 
large number of lower resonances trap disk and halo particles and this trapping is robust, in a broad 
agreement with the earlier results in the literature. 

Subject headings: galaxies: bulges — galaxies: evolution — galaxies: formation — galaxies: halos — 
galaxies: kinematics and dynamics — galaxies: spiral 



1. INTRODUCTION 

Observations and numerical modeling of galactic stel- 
lar bars have been frequently accompanied with basic 
controversies about their origin and evolution. While 
modern understanding of bar growth in a live and re- 
sponsive environment is rooted in the angular momen- 
tum redistribution between the inner and outer disks, 
bulges and dark matter halos (Athanassoula 2003), the 
efficiency of this process is hardly known and its details 
are still to be investigated. Recent efforts include but are 
not limited to the issues related to the bar lifetime cy- 
cles, gas-star interactions, bar amplitudes and sizes, and 
bar slowdown (e.g., Bournaud & Combes 2002; Valen- 
zuela & Klypin 2003; Shen & Sellwood 2004; Weinberg 
1985; Hernquist & Weinberg 1992; Debattista & Sell- 
wood 1998, 2000). The bigger issue of course is how 
the observational and theoretical aspects of bar evolu- 
tion fit within the emerging understanding of cosmolog- 
ical galaxy evolution (e.g., Jogee et al. 2004; Elmegreen 
et al. 2004), specifically the bar evolution in triaxial ha- 
los (El-Zant & Shlosman 2002; Bcrcntzcn, Shlosman & 
Jogee 2005). 

On one hand, early self-consistent models of numer- 
ical stellar bars have relied heavily on the Ostriker & 
Peebles (1973) result which emphasized the Maclauren 
sequence parameter r/|W| — the ratio of bulk kinetic- 
to-gravitational energy as the threshold of bar instability. 
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On the other hand, they revealed robustness of the bars 
— once formed, the bars persisted (Athanassoula 1984 
and refs. therein). The subsequent increase in the par- 
ticle number above TV ~ 10^, switching from 2-D to 3-D 
models with responsive spheroidal components, and in- 
troduction of nonlinear physics tools in the orbital anal- 
ysis have shown a much more complex bar evolution and 
morphology than was anticipated originally (e.g., review 
by Athanassoula 2002a). This refers especially to the 
numerical confirmation that live halos can indeed drive 
the bar instability rather than damp it (Athanassoula & 
Misiriotis 2002; Athanassoula 2003). Lastly, it is still un- 
clear to what extent and how closely the numerical bars 
correspond to their observed counterparts. This issue ex- 
acerbates galaxy studies because overall, both theoreti- 
cally and obscrvationally, the bars appear to be among 
the most important drivers of galactic evolution across a 
wide range of spatial scales. 

In this paper, we have revisited some aspects of a self- 
consistent evolution of stellar bars originating in live stel- 
lar disks embedded in live dark matter halos by focus- 
ing on dynamical and secular changes'* in these systems. 
The evolution of numerical collisionless bars has been 
characterized so far in the literature by three distinct 
phases — the initial growth, the rapid vertical buckling 
and the prolonged quasi-steady regime, i.e., when the 
bars preserve their basic parameters (e.g., Sellwood & 
Wilkinson 1993). However, some indication that bars 
can grow even in the last phase has been noticed already 

We refer to a dynamical evolution when changes develop on 
the timescale of ~ disk rotation, and to a secular evolution when 
they develop on a much longer timescale of ~ 10 — 100 rotations. 
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in low-resolution 3-D models with live halos (e.g., Sell- 
wood 1980). This ability of the bars to grow over ex- 
tended period of time due to the momentum exchange 
with the outer disk and especially with the halo has been 
confirmed recently and analyzed in self-consistent 3-D 
simulations (Athanassoula 2003). Here we attempt to 
quantify this secular growth in terms of the bar size and 
its ellipticity, of the angular momentum exchange, and of 
the ratio of vertical-to-radial dispersion velocities. More- 
over, we look into the corollaries of such recurrent growth 
and find that it leads to additional and substantial 2-D 
and 3-D structural changes in the bar. We, therefore, 
discuss the observational consequences of this evolution. 

Early in their growth stage, numerical stellar bars ex- 
perience a dynamical instability — the vertical buckling. 
The bars thicken profoundly, become more centrally- 
concentrated and acquire a characteristic peanut/boxy 
shape when seen edge-on (Combes et al. 1990; Pfenniger 
& Friedh 1991; Raha et al. 1991; Berentzen et al. 1998; 
Patsis, Skokos & Athanassoula 2002a), while nearly dis- 
solving the outer half of the bar, beyond the vertical in- 
ner Lindblad resonance (Martinez- Valpuesta & Shlosman 
2004). This happens in live models with both axisym- 
metric and mildly triaxial halos (Berentzen, Shlosman & 
Jogee 2005). These boxy/peanut shapes are similar to 
bulge shapes observed in edge-on galaxies (e.g., Jarvis 
1986; Shaw 1987; Bureau & Freeman 1999; Merrifield & 
Kuijken 1999), which can be found in nearly half of all 
edge-on disk galaxies (Liitticke, Dettmar & Pohlen 2000). 
Although observed in numerical simulations a long time 
ago (Combes & Sanders 1981), the origin of boxy/peanut 
bulge shapes still has two alternative explanations — the 
well-known firehose instability (e.g., Toomre 1966; Raha 
et al. 1991; Merritt & Sellwood 1994) and the resonance 
heating (e.g.. Combes et al. 1990; Pfenniger & Friedli 
1991; Patsis et al. 2002b). These two views can be rec- 
onciled if buckling is responsible for shortening the sec- 
ular timescale of particle diffusion out of the disk plane 
and for accelerating the buildup of boxy/peanut bulges 
which proceeds on a much shorter dynamical timescale 
instead (Martinez- Valpuesta & Shlosman 2004). 

However, is the buckling really necessary for a buildup 
of these boxy/peanut shaped bulges? After all, even im- 
posing vertical symmetry did not eliminate this effect, 
albeit the buildup proceeded on a much longer timescale 
(Friedli & Pfenniger 1990). Where are the observational 
counterparts of these asymmetric buckled bars? Due to a 
particular importance of the buckling instability for the 
evolution of numerical bars and its plausible connection 
to the buildup of the pronounced 3-D structure there, 
we have analyzed the bar behavior during and following 
this instability. Specifically, we find that the bars in a 
live environment are capable of recurrent growth, that 
the buckling instability is a recurrent event and that the 
buildup of the 3-D shape is not necessarily a dynamic 
phenomenon. 

Much of the analysis of bar evolution and the accompa- 
ning instabilities is implemented here by means of non- 
linear orbit analysis, because for such strong departures 
from axial symmetry the epicyclic approximation can- 
not be relied upon — being wrong quantitatively it fre- 
quently leads to qualitative errors. We shall try to avoid 
the specific jargon associated with this technique where 
possible. The angular momentum redistribution in the 



model is quantified using orbital spectral analysis. 

In section 2 we provide the details of numerical model- 
ing and analysis. Section 3 describes the overall results of 
secular bar evolution. The 3-D bar orbital structure and 
the inferred vertical structure in the bar are analyzed in 
section 4, and the resonant interaction between the disk 
and the halo in section 5. Discussion and conclusions are 
given in sections 6 and 7. 

2. NUMERICAL TOOLS AND MODELING 
2.1. N-Body Simulations 

To simulate the stellar disk embedded in a live dark 
matter halo, we have used version FTM-4.4 of the iV- 
body code (Heller & Shlosman 1994; Heller 1995) with 
N = 10'"' — 1.1 X 10^. The gravitational forces are 
computed using Dehnen's (2002) f alcON force solver, a 
tree code with mutual cell-cell interactions and complex- 
ity 0(N). It conserves momentum exactly and is about 
10 times faster than an optimally coded Barnes & Hut 
(1986) tree code. 

The initial density distribution is derived from the Fall 
& Efstathiou (1980) disk-halo analytical model. The sys- 
tem is not in exact virial equilibrium, and therefore must 
be relaxed iteratively. The halo-to-disk mass ratio within 
10 kpc is fixed to unity. The halo has a flat density core 
of 2 kpc to avoid excessive stochastic behavior associ- 
ated with the central cusps (El-Zant & Shlosman 2002). 
The disk is exponential and its radial and vertical scale- 
lengths are taken as 2.85 kpc and 0.5 kpc, respectively. 
The disk and halo cut off radii are 25 kpc and 30 kpc 
and the initial circular velocity curves for disk and halo 
components and their sum are given in Fig. 1. The gravi- 
tational softening used is 160 pc and Toomre's parameter 
Q = 1.5. The adopted units are those of G = 1, mass 
M = 10^^ M0 and distance r = 10 kpc. This leads to the 
time unit of T^y^ = 4.7 X 10 yrs and a velocity unit of 
208 km s^^. The energy and angular momentum in the 
system are conserved to within approximately 1% and 
0.05% accuracy, respectively. Above N ~ 10^ our results 
appear to be reasonably independent of N. The model 
evolution presented here has TV = 1.1 x 10^, with 8 x 10^ 
particles in the disk. 
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Fig. 1. — Initial circular rotation velocities for the disk (dashed 
lane) and halo (dotted lane) components. The total is given by the 
solid lane. 

2.2. Orbital and Spectral Analysis 

The self-consistent evolution of stellar bars can be only 
understood by studying their 3-D structure (e.g., Pfen- 
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niger & Friedli 1991; Skokos, Patsis & Athanassoula 
2002a, b) . For such in-depth investigation we use the up- 
dated algorithm described in Heller & Shlosman (1996), 
which is based on a comprehensive search for periodic 
orbits in arbitrary gravitational potentials, and display 
them in characteristic diagrams (section 4) . These orbits 
close in the bar frame and provide the backbone for any 
meaningful analysis of the bar properties. Periodic orbits 
characterize the overall orbital structure of the bar phase 
space, because each of them traps a region of phase space 
around it. These trapped orbits have shapes similar to 
the shapes of the parent periodic orbit. The algorithm is 
run using the package Parallel Virtual Machine (PVM), 
which distributes the search for orbits among different 
processors and computes the orbits using adaptive step 
size, with a relative accuracy of 10^^. Wc also calculate 
the stability of these periodic orbits to estimate which 
orbits can be populated. We track both the 2-D orbital 
families and the 3-D families bifurcating from planar (i.e., 
equatorial) orbits at the (vertically) unstable gaps. The 
shapes of these 3-D orbital families will contribute to the 
shape of the simulated bar, when populated, including 
the evolving shape of the growing boxy/peanut bulge. 

To quantify the disk-halo interaction and the angular 
momentum redsitribution in the system, we have devel- 
oped a package based on the orbital analysis algorithm 
and which uses the Fast Fourier Transform (FFT) to find 
the main orbital frequencies — the angular frequency f2, 
and the radial and vertical epicyclic frequencies k and v 
in the disk and the halo. The particle distribution with 
the frequency ratio is then determined to find the popula- 
tion of resonant orbits. Lastly the change in the angular 
momentum is computed for each of these particles (see 
section 5). 

3. RESULTS 

We first describe the overall bar evolution during the 
simulation period of t ~ 14 Gyr and analyze some of its 
more important behaviors. The bar develops in otherwise 
axisymmetric model during the first two-three disk rota- 
tions. Initially it has an axial ratio (fiatness) c/a ~ 0.1, 
the same thickness as the disk, but this ratio increases 
with time by a factor of ~ 2. The bar starts to brake 
against the outer disk and the halo as seen in Fig. 2b, 
reaching maximum strength at t ~ 1.4 Gyr. At around 
r ~ 1.8 Gyr the bar experiences a vertical buckling in- 
stability which affects its 2-D and 3-D appearance. The 
outer part of the bar nearly dissolves, while overall the 
bar is weakened dramatically, as shown by the to = 2 am- 
plitudes, A2, in Fig. 2a. Immediately following this buck- 
ling, the bar resumes its growth which saturates again at 
r ^ 6 Gyr. At this time A2 in the outer part weakens 
again. The growth is resumed after r ~ 7.5 — 8 Gyr. 
A close inspection of the edge-on bar frames during time 
intervals of 1.8 — 2.8 Gyr and 6 — 7.5 Gyr and the analysis 
described in the next sections, reveal sufficient similari- 
ties between these two events — both represent the ver- 
tical buckling instability in the bar, i.e., the breaking of 
symmetry of the bar (Fig. 3). Such a recurrent buckling 
of bars has never been reported in the literature. 

3.1. Secular Growth of Stellar Bars 

The stellar bar evolution presented here is character- 
ized by substantial changes in the bar size and strength. 
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Fig. 2. — (a). Evolution of bar m = 2 amplitude, A2, for 
r = — 11 kpc (thick solid) and r = 7 — 11 kpc (thin solid) — 
inner and outer bar parts; (b). Bar pattern speed 



and by changes in its 3-D shape (Fig. 3 and Animation 
Sequence 1). Determination of the live bar size in nu- 
merical simulations is not trivial. For example, Athanas- 
soula & Misiriotis (2002) and O'Neil & Dubinski (2003) 
used a variety of methods and found that some of them 
give erroneous and unreliable results. The position of 
A2{r) maximum does not provide any meaningful esti- 
mate for the bar size because the contribution of higher 
harmonics, such as to = 4 and 8, is neglected. To quan- 
tify the bar size changes we have used two alternative 
methods — an ellipse fitting to the isodensity curves and 
the characteristic orbital diagrams. The former method 
has been previously used to detect and to characterize 
'observational' bars (e.g., Knapcn et al. 2000; Laine et 
al. 2002; Hunt & Malkan 2004). Its main deficiency 
when applied to numerical bars is an excessive noise for 
low-to-moderate N and the resulting distribution of the 
bar ellipticities, e(r), being flat for young unbuckled bars 
(Martinez- Valpucsta & Shlosman 2004). However, at 
later times, with the growth of the central mass con- 
centration, this method becomes more reliable. We find 
that a consistently reliable estimate of the bar size, rbar, 
is the radius where e(r) declines ~ 15% from its maximal 
value. The alternative and a new method used by us here 
relates the bar extent to the size of the maximal stable 
orbit of the main orbital family supporting the bar (more 
in Section 4). Both methods produce consistent results 
with each other. 

The bar size evolution is shown in Fig. 4a, for both 
methods. It grows initially to 11 kpc, then buckles and 
shortens to 6 kpc. It then grows again to 13 kpc, where 
the growth stagnates due to the secondary buckling, for 
about 3 Gyr. After this the bar resumes its growth to 
about 16 kpc. The size evolution in the xy plane is ac- 
companied by the vertical thickening of the bar. It does 
not stop after the first buckling instability, but continues 
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Fig. 3. — Evolution of the vertical structure in the bar: edge-on view along the bar minor axis (see also Animation Sequence 1). The 
length is given in kpc and the values of the projected isodensity contours are kept unchanged in all frames. The time in Gyrs is given in 
the upper right corners. The maximal vortical asymmetries correspond to two recurrent bucklings, at t ~ 2.4 Gyr and at ~ 7 Gyr. Note, 
the bar flip-flop between r = 2.3 Gyr and 2.4 Gyr; the persistent vertical asymmetry in the bar at t = 5.2 — 7.5 Gyr; and the development 
of narrow features in the bar midplane outside its core of ~ 8 kpc, after r ~ 9.4 Gyr. Those correspond to ansae ('handles') in the face-on 
bar (see Fig. 11) and are observed in early type disk galaxies. 



gradually due to the vertical resonance scattering, am- 
plified by the recurrent buckling instability of the bar. 
The accompanied boxy /peanut shaped bulge also grows 
with time (Fig. 3). 

With the bar length, rbar, and its corotation radius, 
^CRi we can quantify some of the dynamical character- 
istics of an evolving bar. The ratio rcYi/fhar is shown 
in Fig. 4b. This ratio determines the shape of the offset 
dust lanes in barred galaxies which delineate shocks in 
the gas flow (Athanassoula 1992). The observed shapes 



constrain the ratio rcR/rbar to 1.2 ± 0.2. The modeled 
ratio typically falls within the required limits except dur- 
ing the first buckling when it is higher, 1.5, as noticed 
already by Martinez- Valpuesta & Shlosman (2004). We 
shall discuss this issue in more detail in Section 6. 

3.2. Recurrent Buckling of Stellar Bars 

The buckling is a 3-D phenomena and is most visi- 
ble in the vertical xz-plane (Fig. 3). To quantify the 
bar asymmetry, we have calculated the vertical m = 1 
mode amplitude Aiz in the xz plane and follow the max- 
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Fig. 4. — (a) Evolution of the bar size (semimajor axis). The solid 
line represents the bar size obtained from the ellipse fitting to the 
isodcnsities. The filled dots correspond to the bar sizes obtained 
from the semimajor axis of the last stable x\ periodic orbit which 
supports the bar, obtained from Fig. 7 — a new method introduced 
here (see Section 4 for details); (b) of bar corotation-to-sizc ratio, 
and o{ (c) the maximal ellipticity of the bar, t = 1 ~ b/a from the 
ellipse fitting. 



imal distortion of the bar during its buckling periods. 
Two maxima are apparent at 2.4 Gyr and ~ 7 Gyr in 
Fig. 5. Note, that the vertical asymmetry given by this 
Figure is building up slower and is weaker for the second 
buckling, i.e., from r ^ 5 Gyr to t ~ 7 Gyr the ampli- 
tude is growing to Ai^ ^ 0.03, compared to Ai^ ^ 0.08 
at r ~ 2.4 Gyr. Fig. 3 shows this evolution in a more 
graphical way — while the first buckling affects mostly 
the central few kpc, especially the bar's midplane, the 
second buckling is most prominent in the outer bar range 
of 5-10 kpc and affects the midplane much less visibly. 

The recurrent buckling can be detected in a number 
of ways, e.g., from A2 in the xy-'pla.Tie (Fig. 2a), as men- 
tioned in section 3.1 and above, and from Aiz in the xz 
plane, which quantifies the breaking of vertical symmetry 
in the bar (Fig. 5). Each of the coefficients emphasizes 
a different property of this instability. We are basically 
looking at the same phenomena at different times — the 
first buckling extends over ~ 1 Gyr and the second one 
over ~ 3 Gyr. The changes in the orbital structure of the 
bar during the bucklings will be analyzed in Section 4. 

We have commented in Section 1 that the buckling in- 
stability is a collective breaking of a vertical symmetry 
in the bar. Toomre (1966) has shown that the coupling 
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Fig. 5. — Evolution of the vertical buckling amplitude in the 
bar, i.e. vertical m = 1 mode in the xz plane, integrated over 
r ~ — 12.5 kpc, —2 kpc< y < 2 kpc and —00 < z < 00 intervals. 



between the vertical and radial degrees of motion is the 
prime driver of this instability and the main outcome of 
it is the equalizing of the velocity dispersion in the xy 
plane with the vertical velocity dispersion.^ This evolu- 
tion is characterized by changes in vertical-to-radial ve- 
locity dispersion ratio / . Toomre estimated the crit- 
ical value for this ratio to lie at '^0.1 for a non- rotating 
plane-parallel slab, Raha et al. (1991) at 0.06-0.3 for a 
3-D stellar disk, and Sotnikova & Rodionov (2005) at 
~ 0.6 for the central regions embedded in the hot ha- 
los. Sellwood (1996) have shown that Toomre's limit is 
violated for many of his stellar models, some remaining 
unstable up to 0.4. We plot this ratio for two areas of the 
bar at two different times. For the first buckling (Fig. 6), 
the velocity dispersions are calculated in the central kpc, 
where the maximum effect is expected (e.g.. Fig. 3). Dur- 
ing this buckling, we observe first an increase in dr (the 
initial growth of the bar), and after ~ 1 Gyr a decrease 
in cTr with a corresponding increase in cr^. When tr^/cr^ 
drops below ^ 0.4 the bar buckles and weakens. The 
buckling ends when crl/cr^ increases to unity. The bar 
thickens and grows and the vertical ILR moves gradu- 
ally out. We expect and observe the maximal vertical 
asymmetry in the outer part of the bar during the sec- 
ond buckling, and therefore calculate c^/cr^ at around 
7 kpc. Again, the gradual increase in and decrease in 
(Ti- drive the ratio up from 0.4 to about unity, similar to 
the first buckling. For both bucklings we observe a very 
similar evolution in terms of the velocity dispersions and 
their ratios. 

So far we have shown that the bar buckles twice during 
the simulations: first time — abruptly, fast and in the 
central part, and the second time — slower, less pro- 
nounced and in a different part of the bar. In both 

^ Normally, the kinetic energy of oscillations about the equatorial 
plane, i.e. along the z-axis, is an adiabatic invariant 
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bucklings we observe the loss of symmetry in the vertical 
plane, the drop in A2 and the equilizing of the vertical 
and radial velocity dispersions. After each buckling the 
bar becomes more symmetric, with the difference that 
after the first buckling the asymmetry starts to build 
anew, and after the second buckling — the asymmetry 
is completely washed out. 

To summarize, three main factors characterize the end 
of the buckling instability: A2 starts to decrease, the 
asymmetry in the rz-plane decreases (given by Aiz), and 
(T^/a^ 1. Applying these conditions to our model, the 
end of the first buckling appears at ~ 2.4 — 2.8 Gyr and 
the end of the second buckling at 8 Gyr. 

Since the buckling leads to a sudden vertical thicken- 
ing of the bar and is therefore characterized by particle 
injection above the disk plane, we can monitor this in- 
stability by following the particle distribution. When 
viewed along the bar minor axis, during the first buck- 
ling the bar bends and develops a boxy /peanut shaped 
bulge, while during the second buckling the bar acquires 
an asymmetric shape which leads to the appearance of 
an X-shaped bulge (Fig. 3). These shapes have a direct 
relationship to the population of orbits trapped by the 
main family of periodic orbits in the bar (Pfenniger & 
Friedh 1991; Section 4 below). 

4. ANALYSIS: BAR ORBITAL STRUCTURE EVOLUTION 

Stellar bars are 3-D objects which exhibit dynamical 
and secular evolutionary trends both in their morphol- 
ogy and internal structure. In this section, we study the 
evolution of the orbital structure in the bar by searching 
for the main 2-D and 3-D families of orbits at various 
snapshots, with a particular emphasis on the recurrent 
buckling periods. A comprehensive search for the orbits 
is beyond the scope of this work. We analyze the main 
parameters of detected orbits and calculate their stabil- 



ity. This section uses a specific terminology developed 
for nonlinear orbital dynamics (e.g., Binney & Tremaine 
1987; ScUwood & Wilkinson 1993). Reader unfamiliar 
with this terminology can skip it and go directly to Sec- 
tion 5. We shall discuss the results of the orbital analysis 
in the context of the bar evolution in Section 6. 

The initial vertical axial ratio of the bar which devel- 
ops in our models is c/a ~ 0.1 (nearly a 2-D object), and 
grows dynamically and secularly to ~ 0.2 over the Hub- 
ble time. We first analyze the bar's orbital structure in 
the equatorial plane, then add the 3-D effects accounting 
for the finite thickness and shape of the bar. In doing so, 
we limit our analysis to the periodic orbits within rcR, 
which are largely responsible for the bar shape, and to 
specific times of bar evolution: t = 1.4 Gyr (just prior to 
first buckling), 2.8 Gyr (after the first buckling), 7.1 Gyr 
(during the second buckling) and 11.8 Gyr (after the sec- 
ond buckling). The orbits are searched in the potential 
symmetrized horizontally with respect to the four quad- 
rants. Fig. 5 confirms that the vertical symmetry is vi- 
olated during the buckling periods. The first buckling 
is fast and there is little meaning, therefore, to calculate 
the orbits without some frame averaging — the stars 
see only a time-averaged potential. Hence, we resort to 
a vertically-symmetrised potential at 2.8 Gyr. On the 
other hand, the second buckling is much more gradual 
and we calculate the orbits in the actual 'raw' poten- 
tial at T = 7.1 Gyr to capture the persisting asymmetry. 
The third snapshot with the bar being symmetric again 
is treated similarly to the first snapshot for simplicity. 

The extent of the orbital families and their stability 
are displayed by means of characteristic diagrams. The 
y, z, or z intercept values with the a; = plane are 
plotted with respect to their Jacobi integral, Ej (e.g., 
Binney & Tremaine 1987). The Jacobi (energy) integral 
of motion is conserved along any given orbit in the ro- 
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tating bar frame. The orbits form curves or families in 
the characteristic diagrams. The actual trajectories will 
not coincide exactly with these periodic orbits but may 
be 'trapped' in their vicinity. The properties of periodic 
orbits and their temporal changes, therefore, reflect the 
bar structural evolution. 

The most important single-periodic orbits in the bar 
midplane are those aligned parallel or normal to the bar 
major axis, the so-called Xi and X2 orbits (e.g., Con- 
topoulos & Papayannopoulos 1980; Binney & Tremaine 
1987). The xi constitutes the main family of orbits which 
support the bar figure. While this family exists always 
within rcR , the X2 appear only if the planar inner Lind- 
blad resonance(s) (ILRs) are present. Typically, numer- 
ical bars form with a pattern speed sufficiently high to 
avoid the ILRs, at least in the beginning, and hence the 
xi family completely dominates the bar midplane, short 
of the corotation region. 

The vertical shape of the bar is determined by the 
projection of the populated 3-D orbits onto the corre- 
sponding planes. We search for the vertically-unstable 
gaps in the xi family. It is at these gaps where the 3- 
D families bifurcate through z and z bifurcations (e.g., 
Pfenniger 1984; Pfenniger & Friedli 1991; Skokos et al. 
2002a, b; Patsis et al. 2002b). These gaps coincide with 
the vertical ILR and other vertical resonances. When an 
orbital family goes from being stable to unstable, in the 
mentioned gaps, we get the 3-D prograde families, i.e., 
orbits with initial conditions (y, z, y, z) = (a, b, 0, 0), with 
a, 6 7^ 0. When it goes from being unstable to stable, we 
get the retrograde families, (y, i) = (c, 0,0,d), with 
c, d ^ 0, where x, y and z are oriented along bar's major, 
minor and vertical axes. 

While the vertical ILR is associated with the 
boxy/peanut bulge shapes, which develop as a result of 
the resonance heating of midplane orbits, Patsis et al. 
(2002b) have found that this effect is not limited to a par- 
ticular resonance and not to the barred galaxies per se, 
but operates equally well in nearly axisymmetric and/or 
ovally distorted disks. Furthermore, the buckling itself is 
not a necessary condition for these shapes to form, but 
rather accelerates the process from secular to a dynami- 
cal timescale. 

The main families of 3-D orbits which are responsible 
for the appearance of the boxy /peanut shaped bulges 
are BAN (prograde to the bar rotation) and ABAN (ret- 
rograde) families found by Pfenniger & Friedli (1991). 
They are called xlvl and xlv2, respectively, and appear 
as the 3-D generalizations of xi and X2 families in the 
nomenclature of Skokos et al. (2002a). In general terms, 
the BAN orbits can be described as (z > 0) and 

{z < 0), and ABAN as "00," when projected onto 
the xz-plane. Both families are to : n : ^ = 2 : 2 : 1 
(i.e., two radial oscillations for two vertical oscillations 
for one turn) in the notation of Sellwood & Wilkinson 
(1993). Their projections onto the xy plane have the xi 
orbit shapes. We trace the planar xi and X2 and the 
3-D BAN/ ABAN families from the bar initial growth 
period. For simplicity, we have divided the simulation 
into a number of characteristic time intervals and dis- 
cuss them separately. 

Before and during the first buckling, r ~ — 2.5 Gyr. 
During this period the bar develops and buckles 



(Figs. 2-6). It is geometrically thin and is not centrally 
concentrated, so the vertical ILR is not present. Neither 
present are the planar ILRs and consequently the X2 
orbits. The stable 3-D orbits lie close to the midplane. 
As we stated above, there is little meaning to calculate 
the orbits in the rapidly varying potential without the 
proper time-averaging. Because the midplane of the bar 
is bent in the xz plane (Fig. 3), so are the xi orbits. 

After the first buckling, r ~ 2.5 — 3.5 Gyr. The bar 
has increased its central mass concentration and acquired 
its boxy/peanut shaped bulge (see Fig. 3). As we stated 
above, there is little meaning to calculate the orbits in the 
vertically unsymmetrized potential at the time of the first 
buckling because of the rapidly varying potential. More- 
over, the bar midplane is bent in xz at this time, and 
so are the xi orbits (Fig. 3). Hence, we have vertically 
symmetrized the potential at this time frame (Fig. 7b, 
left frame). In the midplane, the xi orbits dominate and 
no X2 family exists (Fig. 6a). Most of the xi orbits are 
stable, except for a small gap, where the vertical ILR 
is located. At r = 2.8 Gyr the gap is at Ej ~ —2.66, 
where the BAN/ ABAN orbits bifurcate (Fig. 7a). The 
BAN/ABAN are traced (Fig. 7b, c) — their origin in 
z = plane moves gradually out with time toward higher 
Jacobi energy. By the end of the first buckling it stabi- 
lizes and subsequently creeps toward higher Ej as the 
bar brakes against the halo. BAN orbits appear stable 
from 300 pc above the midplane, with a broad unstable 
gap at z ~ 1.1 — 1.3 kpc and a narrow stability island 
(Fig. 7b). Initially, the stable part of these families is 
limited to a very small height above the bar midplane, 
compared to the bar vertical thickness. Gradually the 
stable extent of the BAN family is increasing both in z 
and in x (Fig. 8). The ABAN family is unstable when 
BAN is stable and vice versa (Fig. 7c), as first noted by 
Pfenniger & Friedh (1991). 

Before/during the second buckling, t ~ 3.5 — 8 Gyr. 
The bar is growing both in size and amplitude and 
exhibits a pronounced boxy profile when viewed along 
its minor axis. It appears symmetric with respect to 
the bar midplane up to r ~ 5 Gyr (Fig. 3), although 
Aiz shows that some small residual asymmetry remains 
between the bucklings (Fig. 5). After r ~ 5 Gyr, the 
vertical symmetry of the bar is broken for an extended 
period of time ~ 3 Gyr (Figs. 3, 5). The maximal dis- 
tortion happens at ~ 7.1 Gyr. No X2 orbits exist during 
this time period and the unstable gap in xi moved to 
somewhat higher energies, ~ —2.4 < j ~ —1-7, due to 
the bar braking (Fig. 7a). Additional unstable gap has 
appeared at —2.68 ^ Ej < —2.65, where the 3:2:1 orbit 
family bifurcates. 

Note that the analysis shown in the r = 7.1 Gyr frame 
of Fig. 7 is performed in the actual gravitational potential 
of the TV-body simulation, i.e., without symmetrization 
with respect to the z = plane. Therefore, we should 
not expect to find all the families existing in the sym- 
metrized potentials. But we should be able at least to 
identify the reason for the vertical asymmetry of the par- 
ticle distribution in the projection of the 3D-orbits onto 
the xz plane from the orbital shapes (Fig. 8b). We find 
the BAN family and also the xlv9 family in the nomen- 
clature of Patsis et al. (2002a). 

The BAN family appears stable everywhere above the 
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Fig. 7. — Three snapshots of the bar orbital structure evolution. Left — immediately after the first buckling, middle — during the 
second buckling, right — well after the second buckling. The t = 7. 1 Gyr frame is non-symmetrized (vertically), while others are vertically 
symmetrized, (a). Characteristic diagrams in the xj/-plane showing j/-intersections of regular orbits with the a; = axis as a function of 
Jacobi energy of the orbits. The bar is oriented along the x-axis. Solid lines represent stable orbits, dotted — unstable ones. The dashed 
lines show the ZVCs — the zero velocity curves, (h). Characteristic diagram in the a:2:-plane showing ^-intersections of regular orbits. 
Since 7.1 Gyr snapshot is non-symmetrized and the bar buckles in this plane, the distribution of the BAN orbits is also asymmetric, and 
most of the stable BAN orbits are (see text). They have a larger curvature and extend further into the 2 < region leading to the 

bar asymmetry in the edge-on view (Fig. 3). (c). Characteristic diagrams in the 2;2:-plane showing the z families. The 7.1 Gyr snapshot is 
omitted because all the ABAN orbits are unstable at this time and, therefore, difficult to calculate in unsymmetrized potential. 



midplane, except around z ~ 2 kpc (Fig. 7b), while 
ABAN is unstable everywhere and so is not shown here. 
For z ^ 0; the BAN family is detached from the mid- 
plane and starts at z ~ —0.5 kpc. It forms a kind of loop 
in the characteristic diagram. We do not find any stable 
orbit below z ~ —2 kpc. The x\v% family is stable only 
for a low range in energies and only below the midplane, 
which makes it important for contributing to the vertical 
asymmetry in the bar (Fig. 8b). 

After the second buckling, t > 8 Gyr. The bar slowly 
regains its vertical symmetry, and becomes symmetric 
after ~ 8 Gyr. As a result of this instability, the bar 
has developed a pronounced X-shape which differs from 
the 'usual' peanut shape because of the extended con- 
cave region between the spikes, i.e., peanuts (Fig. 8c). 
Orbits which appear at the same z in the characteristic 
diagram but at different time frames differ in their pro- 
jections. After the second buckling the orbits are longer 
and extend higher, i.e., the convexity and concavity have 
significantly changed. The vertical ILR has moved to 

Ej 2.36 by 11.8 Gyr. The BAN orbits are stable 

while ABAN are unstable. Overall, the boxy/peanut 
bulge has increased in size from ~ 1 — 2 kpc (after the 
first buckling) to ~ 8 kpc (after the second buckling). 

The X2 orbits and consequently two ILRs in the xy- 



plane appear only after the second buckling, > 7.2 Gyr. 
Initially occupying a small range in energies, their ex- 
tent grows with time. The xi family is mostly stable up 
to the (mostly) unstable 'shoulder' in the characteristic 
diagram (Fig. 7a, right frame), which has continued to 
move out to higher Jacobi energies, along with the last 
stable orbit supporting the bar. 

We also use the characteristic diagrams to get an in- 
dependent estimate for the bar physical size. Unlike the 
case of analytical bars where this issue is resolved triv- 
ially, to estimate the length of the live numerical and 
'observed' bars can be more difficult. For example, if the 
bar potential or bar nonaxisymmetric force is used for 
this purpose, one can get an erroneous result that the 
bar extends beyond its corotation radius. If one relies on 
the density distribution, one can overestimate the size as 
well because the bar can drive a pair of open spirals. In- 
stead, we use the calculated orbital structure of the bar 
and rely on the properties of Xi family which is generally 
stable except for narrow gaps (in Jacobi energy) and in 
the region around the corotation. Specifically, close to 
the corotation, the xi curve bends upward (e.g.. Fig. 7a; 
Binney & Tremaine 1987). We use the x-axis extent of 
the last stable xi orbit which lies on this upward branch. 
It is the slow drift of unstable gaps and of the elbow of the 
si-curve towards higher E] which leads to the trapping 
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Fig. 8. — Projection of 3-D orbits onto the xz plane in three 
different snapshots — after the first buckling, during and after the 
second buckling. Only stable orbits of the main families described 
in the text and in Fig. 7 arc plotted. The horizontal extension 
of the orbits increases with time. The general shape evolves from 
boxy/peanut to a vertically-asymmetric, and then to a peanut/X 
shape. 



of orbits and the secular increase in the bar size. A large 
number of models here and in Berentzen et al. (2005) 
have been used for the comparison with the isodensities 
fit to the bar (Section 3). 

The bar amplitude A2 shows that the outer half of the 
bar practically dissolves in the first buckling and weakens 
substantially in the recurrent buckling (Fig. 2). The el- 
lipse fit to the bar reflects this by a dramatic decrease in 
the bar length from 11 kpc to ~ 6 kpc in radius. Un- 
fortunately, it is more tricky to estimate the bar size from 
the orbital analysis at this time — the underlying poten- 
tial is time-dependent while the analysis is performed in 
the frozen potential. Nevertheless, both methods agree 
before and after the the first buckling and capture the 
change in the growth rate of the bar length after the 
second buckling. To summarize, the bar size decreases 
during the first buckling and levels-off during the second 
one. 

The simple explanation for the success of the orbital 
analysis method is that it is the most self-consistent 
method known to us and relies on the orbit trapping by 
the stable periodic family whose existence basically de- 
fines the bar. We note that the middle frame of Fig. 7a 
does show a large unstable gap just below the elbow of xi 
curve. However, we note once more that the potential in 
this frame has not been vertically symmetrized and that 



the weakening of the bar during this time of buckling is 
in fact a direct consequence of this gap. 

In summary, the characteristic diagrams and their re- 
spective orbital families reflect both the dynamical and 
secular evolution of the bar shown in both Fig. 3 and 
the Animation Sequence. The ustable gaps, where the 
vertical families (BAN/ABAN) bifurcate, move toward 
higher Jacobi energies while the orbits become more ex- 
tended along the bar axes. The orbital families evolve not 
just quantitatively but also qualitatively — this shows 
up in the changing shape of the bar. Although we do 
not measure the population of the different families, the 
bar shapes appear to be governed mainly by the BAN 
family, which evolves by changing the concavity of its or- 
bits, thus giving the bar the pronounced peanut/boxy/X 
shapes when viewed along the minor axis. The change 
in the orbital structure of the bar supports our obser- 
vation, based on independent arguments, that there is a 
recurrent buckling. At the secondary buckling time, the 
vertical families differ above and below the plane. The 
former are stable orbits extending more along the x-axis 
and also giving it more of a boxy shape. The latter ones 
are stable orbits which give the bar more of a peanut 
shape. The appearance of the X2 family confirms that 
the ILRs form only after the second buckling of the bar. 

5. DISK-HALO ANGULAR MOMENTUM EXCHANGE: THE 
ROLE OF THE RESONANCES 

The disk region which gives rise to a bar is known to 
lose its angular momentum (J). In principle, a number 
of components in a galaxy can acquire this momentum, 
namely, the outer disk, the bulge, and the halo. While 
the former two components are able to store J, the re- 
sponsive (i.e., live) halo can serve as a particularly large 
angular momentum sink due to its large mass and low J. 
Athanassoula (2002b, 2003) has shown that the disk-halo 
interaction is mediated by lower resonances. In order to 
understand more fully the complex behavior exhibited by 
the model, such as the initial bar growth, its recurrent 
bucklings, and subsequent secular growth — we find it 
instructive to determine the balance of the angular mo- 
mentum and its evolution. We pay particular attention 
to the resonant interactions between the model compo- 
nents. 

We first determine the overall J balance in the different 
regions of the disk and the halo. Fig. 9 shows the total 
angular momentum evolution for three different regions 
in the disk and halo: inner (0 — 7 kpc), intermediate (7 — 
15 kpc) and outer (15 — 35 kpc or 15-50 kpc) regions. The 
inner and intermediate regions are those hosting the bar 
for various times during its evolution. The time of the bar 
formation is characterized by a substantial loss of J in the 
inner disk. This angular momentum is redistributed to 
the intermediate and especially to the outer disk (Fig. 8). 
During the first buckling, there is some indication that 
the inner disk gains some angular momentum, and we 
return to this issue in section 6.2. The subsequent growth 
of the bar is accompanied by a slow J growth in the outer 
disk and the halo. Interestingly, after the second buckling 
J stops growing in the outer disk, while the outer halo 
picks-up J at an increasing rate. There is a general loss of 
angular momentum from the disk and there is a general 
gain in the halo. The loss of J in the disk correlates with 
the slowdown of the bar. which is known to anti-correlate 
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Fig. 9. — Evolution of angular momentum in the various regions 
of the disk (upper) and the halo (lower). The specified radii are 
given in cylindrical geometry. 



with the growth of the bar (e.g., Athanassoula 2003). 

Next, we calculate the contribution of the resonances 
between disk and halo particles with the bar's f2b to the 
angular momentum exchange between the disk and the 
halo. We generally follow the procedure described in 
Athanassoula (2002b). The principal frequences, fl and 
K of azimuthal motion and of radial oscillations of the 
3-D orbits, have been determined by means of spectral 
analysis (Binney & Spergel 1982). We freeze the gravi- 
tational potential in the model at times t = 2.8 Gyr and 
5.2 Gyr, but allow the bar to tumble with its original Q]^, 
and integrate the orbits of 1.5 x 10'"' randomly-picked par- 
ticles in the disk and the same amount in the halo for an 
additional 9.4 Gyr, about 20 bar tumblings. Results are 
plotted in Fig. 10 as a function of the frequency ratio 
?7 = (f2 — rib)/^, where TVr particles are binned in the 
intervals of Ar; ~ 0.005 (top frames). The difference A J 
in angular momentum between the two times is given as 
a function of ry at r = 5.2 Gyr (bottom frames). The ma- 
jor resonances are located at ?/ = ±0.5 (ILR and OLR), 
±1 : 3, 1 : 4, etc., with the positive rj indicating the reso- 
nances inside the corotation radius and ?7 < — outside 
the corotation. The corotation resonance corresponds to 
?7 = 0. 

The resonant and near-resonant particle distributions 



N-n at the time t = 2.8 Gyr are clearly non-uniform and 
are permeated by numerous resonances. The trapping 
of orbits by the resonances is also evident. The domi- 
nant resonance in the disk appears to be the ILR which 
traps the largest number of particles and therefore facil- 
itates the loss of the angular momentum from the disk 
(lower panels of Fig. 10). Note, that we have cut 
for this resonance to keep the same scales for disk and 
halo diagrams — it extends to 39 x 10'^. The xi within 
the CR and later on X2 (if populated) orbits, as well as 
BAN/ABAN orbits are trapped within this resonance. 
The broad asymmetric peak (Fig. 10a, left panel) in the 
range of 77 ~ ±0.35 — 045 is made out of near-resonant or- 
bits in the (rotating) boxy/peanut bulge, many of them 
near resonant BAN orbits. Those have been injected 
during the first buckling. They appear on both sides 
of the CR because J along these orbits oscillates wildly, 
A J ~ J, and it is sometimes difficult to disentangle the 
prograde from retrograde ones among them. For com- 
parison, we repeated this procedure at r = 11.8 Gyr 
and confirmed that more than half of the particles in 
the above broad peaks have been trapped by the nearby 
ILR and OLR by that time. The resonances for 77 ^ 0, 
with CR and OLR as the next strongest ones, absorb the 
angular momentum and their A J are all positive. 

The halo is dominated by the corotation resonance 
which absorbs J from the disk. Overall, the disk reso- 
nances emit and the halo resonances absorb J. This be- 
havior was demonstrated by Athanassoula (2002b, 2003), 
where similar plots to Fig. 10 were given and there is a 
very good agreement between our results. The only dif- 
ference we find is that the halo ILR is actually losing J 
in our model, albeit a small amount, unlike in Athanas- 
soula's models. The modeled halo extends well beyond 
the disk and its outer part appears to be actively storing 
J, up to ?/ ^ — 3. This effect is confirmed by the overall 
J evolution in the halo (Fig. 9). We find that the trap- 
ping of orbits by the resonances is robust — almost all 
the particles trapped at t = 2.8 Gyr remain trapped at 
5.2 Gyr. The ability of the barred disk to transfer its 
angular momentum to the halo while J saturates in the 
outer disk, explains why the bar is able to grow during 
this time interval as exhibited by the A2 amplitude and 
the size of the bar (e.g.. Figs. 2-5). It is this growth of 
the bar which is ultimately responsible for the recurrent 
buckling phenomenon analyzed here. 

6. DISCUSSION 

We have modeled the evolution of coUisionless (stellar) 
bars embedded in responsive (live) axisymmetric dark 
matter halos using high-resolution iV-body simulations. 
We find that (1) bars experience secular growth over 
the simulation (~ Hubble) time, except (2) during well 
defined time periods when the bars encounter sponta- 
neous breaks of vertical symmetry, so-called buckling (or 
firehose) instability (e.g.. Fig. 3 and the Animation Se- 
quence 1). We detect such a recurrent buckling during 
which the bar, and especially its outer half, weakens sub- 
stantially, but the growth is resumed subsequently. Two 
different techniques have been used to measure the bar 
size, namely, the nonlinear orbital analysis and the iso- 
dcnsity ellipse fitting. The bar strength has been mea- 
sured using the amplitude of m = 2 mode, and the bar 
vertical asymmetry — using the vertical m = 1 mode. 
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Fig. 10. — Resonance interaction in the live disk-halo system. Upper: histogram of particle distribution N-n with the principal frequency 
ratio ri = {Q — Qi,)/k for the disk and (rotating) boxy/peanut bulge (left) and halo (right), about 1.5 X 10^ randomly picked particles each, 
at r = 2.8 Gyr. The cusp at the ILR (r] = 0.5) for the disk is shown only partially — it extends to ~ 39 X 10^. The broad peaks near 
r] ~ ±0.4 are made mostly of near-resonant bulge particles injected during the first buckhng. Lower: angular momentum difference, A J, 
between the particles at t = 5.2 Gyr and t = 2.8 Gyr in the disk (left) and halo (right) as a function of rj at time 5.2 Gyr. The lower 
resonances, such as ILR [rj = 0.5), the corotation {rj = 0), OLR (rj = —0.5), 1:3 (/y = ±0.33), 1:4 (rj = 0.25), etc. are clearly visible. The 
ILR dominates the loss of the angular momentum in the disk, and the corotation dominates the absorption of the momentum by the halo 
(more in the text). The frequency ratio is binned into intervals of Ar; = 0.005. Note that r] scales are different for the disk and the halo. 



Moreover, we have analyzed the bar evolution and its 
buckling periods by means of the nonlinear orbit analy- 
sis and the ratio of vertical-to-planar velocity dispersions. 
Finally, we have examined the bar development in terms 
of the angular momentum redistribution between vari- 
ous components in the disk-halo system mediated by the 
resonant interactions using orbital spectral analysis. 

Low resolution iV-body simulations exhibit fast (e.g., 
during one disk rotation) stellar bar growth, followed by a 
vertical buckling and secular weakening. High-resolution 
simulations with N 10^^^ show a more complex evo- 
lution. They allow for modeling the disks embedded in 
live halos and account for the resonant interactions be- 
tween the halo and the disk, adding a new and crucial 
element to the simulations in the form of angular mo- 
mentum transfer between these components. The main 
difference with the previous low-resolution models and 
that presented in this work is the ability of the stellar 
bar to strengthen again after its original weakening fol- 
lowing the buckling — a process which leads to the bar 
growth and consequently to its recurrent buckling. 



6.1. Disk-Halo Resonant Interactions 

The bar growth has been associated with the exis- 
tence of 'sinks' of angular momentum located elsewhere. 
Athanassoula (2002b, 2003) has shown that a galactic 
halo can play such a role and absorb large quantities of 
angular momentum from the disk/bar region. This effect 
appears to completely invert the original suggestion by 
□striker & Peebles (1973) about the stabilizing function 
of dark matter halos in disk galaxies against the bar for- 
mation instability. Angular momentum redistribution, 
at least in principle, can have contributions both from 
the non-resonant and resonant interactions between the 
bar and orbits in the halo. To capture the latter requires 
large N to stabilize the population of resonant particles 
(e.g., Weinberg & Katz 2002). 

It is therefore important to demonstrate that resonant 
particles are indeed responsible for the angular momen- 
tum transfer in the first place. We start with the relative 
contributions to the J transfer between the disk and the 
halo and between the inner (i.e., bar unstable ~ 15 kpc) 
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and outer disks (Fig. 9). Although the efficiency of this 
redistribution is of course model dependent, we never- 
theless quote the numbers, assuming that to a certain 
degree, they are representative. We find that ^ 39% of 
the angular momentum in the disk is lost to the halo 
during the evolution. The bar unstable region, <, 15 kpc 
has lost about 71% of its original J, of which about 38% 
went to the outer disk and the rest was absorbed by the 
halo. The most intensive flow of J to the outer disk hap- 
pens in the early stage which ends with the first buckling. 
The halo particles have much larger dispersion velocities 
than the disk and until the bar fully develops, their res- 
onant interaction with the m = 2 asymmetry in the disk 
is virtually non-existent. 

After the first buckling, it is the halo which absorbs 
most of the angular momentum lost by the bar region 
and from Fig. 9 (lower frames) it is clear that the lion 
share of this exchange is resonant and mediated by the 
CR and the OLR in the halo and by the ILK and the CR 
in the disk, with some contribution from the numerous 
minor resonances. This results in the increase of the bar 
size and strength, and leads to the gradual decrease in the 
ratio of the velocity dispersion, /cr^ , in a close analogy 
with the disk evolution preceeding the first buckling. The 
secondary buckling of the bar, therefore, can be directly 
traced to its robust growth in our model. 

The issue of a bar dissolution vs. growth has become 
controversial recently. While it is beyond the scope of 
this work, indirectly we do touch it. First, we confirm the 
results of Martinez- Valpuesta & Shlosman (2004) that 
vertical bucklings do not destroy the bar, unlike sug- 
gested by Raha et al. (1991; see also Sellwood & Wilkin- 
son 1993). Next, in various axisymmetric models, here 
and in Berentzen et al. (2005), a robust growth in the 
bar is observed after the first buckling, supplemented by 
the angular momentum exchange between the disk and 
the halo. Athanassoula (2003) has analyzed the effect 
of initial conditions on this J redistribution, namely, of 
the initial disk-to-halo mass ratio and of the Toomre's 
Q-parameter — the angular momentum exchange varied 
from ~ 2% to ~ 35% with a clear correlation between 
the amount of the exchanged momentum and the size of 
the bar. In comparison, our models have a halo-to-disk 
mass ratio of unity within the central 10 kpc, Q ^ 1.5, 
and a halo core of ^ 2 kpc. They lie between the models 
MHl, M73-4 and MQ4 of Athanassoula. 

Interestingly, Valenzuela & Klypin (2003) argued that 
it is the mass and force resolutions of numerical models 
which dictate the efficiency of J transfer to the halo and 
hence play an important role in the evolution of the bar. 
Insufficient resolution results in an excessive growth of 
numerical bars and an excessive decrease in the bar pat- 
tern speed. However, it seems rather that the physical 
conditions mentioned above take priority — their high 
resolution models have extended and hot halos. High dis- 
persion velocities in the halo affect the particle trapping 
by the resonances and so are expected to lower substan- 
tially the J transfer. 

6.2. Bar size evolution 

If one neglects the dissipative component in a modeled 
galaxy, the bar pattern speed decreases secularly, except 
during brief time intervals of internal instabilities in the 
bar itself. At least in numerical models of stellar bars. 



their overall slowdown is accompanied by an increase in 
the bar length, so the bar roughly extends to its coro- 
tation radius which increases with time (Athanassoula 
1992). Alternative theoretical models of bars terminat- 
ing at the ILR exist (Lynden-Bell 1979), but difficulties 
remain in actually reproducing them in numerical sim- 
ulations (but see e.g., Polyachenko & Polyacnenko 1994 
for the so-called slow bars). A sole exception consists of 
a system of nested bars, where the inner (nuclear) bars 
form within the ILR (e.g., Shlosman, Frank & Begelman 
1989; Friedfi & Martinet 1993; Englmaier & Shlosman 
2004). However, a dissipative component is required to 
be present for self-consistency in this case. How exactly 
the bar traps the disk orbits in order to increase its length 
is not known at present. 

We have demosntrated that while the size of the bar, 
'"bar, drops during the first buckling (Fig. 2a), the ra- 
tio of rcn/rbar (Fig. lb) does not increase dramatically 
above the range of 1.2 ±0.2 determined by Athanassoula 
(1992) to fit the observed shapes of the offset dust lanes 
in barred galaxies. Amazingly, the reason for this is that 
rcR drops as well due to the sudden increase in the bar 
pattern speed in the same time interval. This is not unex- 
pected due to a clear trend between fib and A2 (Athanas- 
soula 2003) and a long known fact that the bar lengthens 
with its slowdown in numerical simulations. What is new 
here is that we find that f2b increases sharply during the 
first buckling and this increase apparently correlates with 
the sudden decrease in the bar length. Hence the rela- 
tion between between fib, A2 and rbar appears to be more 
fundamental than anticipated and holds for decelerating 
and accelerating bars altogether. 

We associate the continued bar growth with the abil- 
ity of the halo to absorb angular momentum from the 
disk region lying within the corotation radius — a region 
which itself expands with time. This process leads to 
a stronger bar and can be seen as a counterbalance to 
a number of other processes which have been discussed 
in the literature within the framework of bar dissolution. 
Recent observational results from Galaxy Evolution from 
Morphology and SEDs (GEMS, Rix et al. 2004) sur- 
vey have shown that the bar size and axial ratio dis- 
tributions at intermediate redshifts of z ~ 0.2 — 1 are 
compatible with those in the local Universe (Jogee et 
al. 2004; Elmegreen et al. 2004; Sheth et al. 2003). 
Our models which exhibit a slow secular bar growth over 
the Hubble time are in general agreement with these re- 
sults. They do lack the gaseous component which can 
dramatically shorten the bar life cycle (e.g., Bournaud 
& Combes 2002). However, it is difficult to understand 
how the short-lived (e.g., less than 1.5 Gyr [Bournaud & 
Combes 2002; Combes 2005]) bars can form at the same 
rate and with the same sizes and axial ratios they are 
being destroyed by the gas. Yet, taken at the face value, 
in order to agree with the GEMS results, the bars must 
preserve their size and strength distributions nearly un- 
changed over the last 8 Gyr. Clearly, the precise role of 
gas in the bar evolution must still be determined. 

6.3. Bar shape evolution 

Next we focus on some aspects of the bar 3-D shape 
evolution in the model, specifically on the axial ratio in 
the bar midplane and on the bar symmetry in the ver- 
tical plane. The midplane axial ratio, or alternatively 
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the bar ellipticity e, is measured from the ellipse fitting 
(Fig. 4c) and agrees with the evolution of the bar ampli- 
tude A2. At the same time, the variations in e appear 
less dramatic than in A2, although recurrent bucklings 
are clearly visible as a decrease or saturation. 

The ratio rcn/rhar has also dynamical implications for 
the bar. For example, it determines the shape of the 
offset dust lanes in barred galaxies, which in turn delin- 
eate shocks in the gas flow (Athanassoula 1992). The 
strength of the underlying shocks determines the gas in- 
flow towards the central kpc. A bar which is considerably 
weakened will slow down the radial gas inflow. The ob- 
served shapes constrain the ratio rcn/rhar to 1.2 ± 0.2. 
The modeled ratio (Fig. 4b) typically falls within the re- 
quired limits except during the first buckling when it is 
higher, ~ 1.5. 
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Fig. 11. — Ansae at time t = 9.9 Gyr. A face-on view of the 
modeled barred disk. The gray-scale contours represent the surface 
density and are spaced logarithmically. The box size is 18 kpc X 
18 kpc. 

After the second buckling, the bar shows the so-called 
ansae (handles) on both its ends. We observe them 
as characteristic density enhancements in the face-on 
(Fig. 11) or edge-on (Fig. 3) disks. Athanassoula (2001) 
related the appearance of ansae to initial conditions in 
the models (e.g., the halo-to-disk mass ratio). It remains 
unknown why they appear at this particular evolution- 
ary stage, T > 8 Gyr. The ansae can be seen in some 
early-type barred galaxies, e.g., NGC 4262, NGC 2859 
and NGC 2950 (Sandage 1961), NGC 4151 (Mundell & 
Shone 1999), ESQ 509-98 (Buta et al. 1998), etc. 

Vertically, the bar evolves from a geometrically thin 
configuration, similar to the disk hosting it. The verti- 
cal bar buckling, when viewed along the bar's minor axis, 
shows a rapidly evolving bending which relaxes to a boxy 
shaped bulge, i.e., bulge with flat or mildly convex iso- 
densities (Fig. 8a). Subsequently, the bulge acquires a 
peanut shape. During the second buckling, for ~ 3 Gyr, 
the vertical asymmetry persists with one-sided boxy and 



peanut symmetries which derive from the asymmetric 
BAN family (Fig. 8b). With the asymmetry washed out, 
the boxy/peanut bulge/bar has a pronounced X-shape: 
two pointed spikes with a large concave region in be- 
tween (Fig. 8c). The X-shapcs have been seen before in 
numerical simulations (e.g., Athanassoula & Misiriotis 
2002) and have been observed as well (e.g., NGC 4845, 
NGC 1381, IC 4767 (Whitmore & Bell 1988), IC 3370 
(Jarvis 1987) and AM 1025-401 (Arp & Madore 1987; 
Bureau & Freeman 1999; Patsis et al. 2002a). Mihos 
et al. (1995) proposed a merging scenario for the for- 
mation of the X-shaped bulges, e.g., a minor merger for 
Hickson 87a galaxy. This merger triggers the bar insta- 
bility in the disk, followed by the buckling and the X- 
shaped bulge, when viewed from a specific aspect angle. 
In principle, there may be more than one cause for the 
bar buckling and for the formation of peanut/boxy/X- 
shaped bulges. But when taken together, a high observed 
frequency of these bulges (Liitticke et al. 2000) and the 
need for only mild asymmetry in the disk for their ap- 
pearance (Patsis et al. 2002b) may hint to their intrinsic 
origin. 

6.4. Vertical asymmetry of the bar 

The vertical asymmetry of the modeled bar has been 
detected during its recurrent bucklings, for about 1 Gyr 
and ~ 3 Gyr, respectively (Figs. 3, 5). Its characteristic 
shape, in principle, is detectable by observations of edge- 
on galaxies (e.g.. Fig. 3 and the Animation Seq. 1), espe- 
cially during the second buckling due to its prolonged pe- 
riod. Observationally, the bucklings differ when viewed 
edge-on along the bar minor axis. The main difference 
is the location of the maximal asymmetry — it is close 
to the rotation axis during the first buckling and around 
the middle region of the bar during the second buckling. 

Liitticke et al. (2000) present first statistics of edge- 
on galaxies with the boxy/peanut bulges and show some 
isophote fits (their Fig. 2). They distinguish between 
peanut-shaped and boxy-shaped bulges. A total of 27% 
of the sample of 734 galaxies with boxy/peanut bulges 
have bulges which are either close to boxy or, due to low 
resolution, could not be distinguished between boxy and 
peanuts. Some of these bulges show a mild vertical asym- 
metry similar to that found during the second buckling, 
e.g., NGC 4289 — a good example of a possible observed 
secondary buckling. Of course, to corroborate the numer- 
ical simulations one needs a large, statistically significant 
sample with high resolution. The main difficulty lies in 
the unambigious determination of the inclination angle 
of a galactic disk — it must be edge-on within ±5° — 7°. 

From a theoretical point of view, it is unclear how 
wide-spread is the recurrent buckling in stellar bars — 
how sensitive it is to disk and especially halo parame- 
ters, such as mass distribution and dispersion velocities. 
While the condition for a second buckling seem to be di- 
rectly related to the ability of the bar to strengthen after 
the initial weakening, this process can depend on a num- 
ber of additional parameters. For example, Athanassoula 
(2003) has found that the angular momentum transfer 
between the disk and the halo depends on the mass dis- 
tribution in the halo and weakens substantially for hotter 
halos. Berentzen et al. (2005) have detected a recurrent 
buckling in their LS2 model with 2 kpc fiat core (live) 
halo with a logarithmic potential. At the time of the 
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second buckling, the bar in this model appears some- 
what stronger than in the first buckling. Moreover, the 
secondary buckling can be apparently seen in Fig. 4 of 
O'Neil & Dubinski (2003) and also in Fig. 11 of Valcn- 
zuela & Klypin (2003), the latter based on the evolution 
of the bar pattern speed, strength and the angular mo- 
mentum rate change. Secondary buckling is also present 
in some models of L. Athanassoula (private communica- 
tion 2005). In all the cases they went unnoticed. Lastly, 
we comment on the addition of the dissipative component 
to the stellar disk (Berentzen et al. 1998). The effect 
of the clumpy gas component is to weaken the buckling 
instability, but without quantifying the degree of dumpi- 
ness within the central kpc, it is not clear whether the 
isothermal equation of state used has led in fact to over- 
damping of this instability. 

7. CONCLUSIONS 

To summarize this work, we have studied the long-term 
stellar bar evolution in a high-resolution self-consistent 
model of a disk and a responsive halo. Wc find that a 
developing bar goes through the vertical buckling insta- 
bility which weakens it and dissolves its outer half. Sub- 
sequently, the bar experiences a renewed growth which 
leads to a recurrent buckling. This evolution is driven by 
the resonant interaction between the barred disk and the 
surrounding halo — we quantify this effect by means of 
the spectral analysis of individual orbits in the disk and 
the halo and show that the halo particles are trapped by 
numerous lower resonances with the bar and that this 
trapping is robust. During these periods of recurrent 
instability, and especially during the slower second buck- 
ling, the bar remains vertically asymmetric for a pro- 
longed ^ 3 Gyr time interval — which in principle can 
be detected observationally. However, two issues can po- 
tentially complicate this detection. First, it is not clear 
how widespread are the conditions favorable for the re- 
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current bar growth, although wc have detected it in a 
number of models with different initial conditions. Sec- 
ond, while it was shown that a clumpy gaseous compo- 
nent with an isothermal equation of state in the disk 
will weaken this instability, the effect of a realistic ISM 
was never estimated. A statistically significant sample of 
(nearly) edge-on galaxies is required to test the predic- 
tion of a prolonged vertical asymmetry. 

Wc also find that the secular bar growth and the trig- 
gered buckling instabilities lead to pronounced changes 
in the bulge shape — it grows both radially and ver- 
tically, acquiring a peanut, a boxy and finally the X- 
shaped appearance. While the bar size approximately 
follows the 4:1 (Ultra-Harmonic) resonance in the disk, 
the boxy/peanut bulge size appears to be guided by the 
vertical ILR. Concurrently, the bar is going through a 
structural evolution — new families of 3-D periodic or- 
bits appear (or become more pronounced) after the buck- 
lings. 

Finally, wc find that the bar strength correlates with its 
pattern speed, in both strengthening and weakening bars. 
While it was known already that a stellar bar becomes 
stronger as it slows down, we detect the reverse trend 
as well — bars that weaken during the buckling speed 
up their tumbling. Moreover, the bar size appears to be 
sensitive to these changes — slowing down bars become 
longer, while speeding up bars shorten. 
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